N_temp = 1000;
ngrid_x = 25;

ratio_temp = NaN(N_temp,1);
increment_x = NaN(N_temp,1);
discr_perc = NaN(N_temp,1);
ratio_temp_first = NaN(N_temp,1);
increment_x_first = NaN(N_temp,1);
discr_perc_first = NaN(N_temp,1);
monotone = NaN(N_temp,1);

for i=1:N_temp

    [~,temp] = max(data.z(i,:));

    temp_diff = setdiff(1:J,temp);

    ind_good2 = temp_diff(1);

    argument = [data.x(i,temp) data.x(i,temp_diff) data.z(i,temp) data.z(i,temp_diff)];
    argument_dz1 = argument;
    argument_dz1(J+1) = argument_dz1(J+1) + increment_z;
    argument_dz2 = argument;
    argument_dz2(J+2) = argument_dz2(J+2) + increment_z;
    argument_dz1_dz2 = argument_dz1;
    argument_dz1_dz2(J+2) = argument_dz1_dz2(J+2) + increment_z;

    s = fun_sigma_est(argument,theta_NPD,combos_all_unique,m_regr);
    s_dz1 = fun_sigma_est(argument_dz1,theta_NPD,combos_all_unique,m_regr);
    s_dz2 = fun_sigma_est(argument_dz2,theta_NPD,combos_all_unique,m_regr);
    s_dz1_dz2 = fun_sigma_est(argument_dz1_dz2,theta_NPD,combos_all_unique,m_regr);

    d2s_dzdz = s_dz1_dz2-s_dz2;

    grid_x = linspace(data.x(i,temp_diff(1)),1,ngrid_x);
    d2s_dzdx_temp = NaN(ngrid_x,1);
    ds_dx_temp = NaN(ngrid_x,1);
    for xx=1:ngrid_x
        argument_dx2 = argument;
        argument_dx2(2) = grid_x(xx);
        argument_dz1_dx2 = argument_dz1;
        argument_dz1_dx2(2) = grid_x(xx);

        s_dx2 = fun_sigma_est(argument_dx2,theta_NPD,combos_all_unique,m_regr);
        s_dz1_dx2 = fun_sigma_est(argument_dz1_dx2,theta_NPD,combos_all_unique,m_regr);

        d2s_dzdx_temp(xx) = s_dz1_dx2-s_dx2;
        ds_dx_temp(xx) = s_dx2-s;
    end

    % if all(diff(d2s_dzdx_temp)>0) || all(diff(d2s_dzdx_temp)<0)  % only do this if d2s_dzdx_temp is monotonic

    [discr,ind_x] = min(abs(d2s_dzdx_temp-d2s_dzdz));
    discr_perc(i) = abs(discr/d2s_dzdz);
    monotone(i) = (all(diff(d2s_dzdx_temp)>0) || all(diff(d2s_dzdx_temp)<0));
    increment_x(i) = grid_x(ind_x)-data.x(i,temp_diff(1));
    ratio_temp(i) = increment_x(i)/increment_z;

    [discr_first,ind_x_first] = min(abs(ds_dx_temp-(s_dz2-s)));
    discr_perc_first(i) = abs(discr_first/(s_dz2-s));
    increment_x_first(i) = grid_x(ind_x_first)-data.x(i,temp_diff(1));
    ratio_temp_first(i) = increment_x_first(i)/increment_z;

    %i
end

display('*** SECOND DERIVATIVES ***')
select = (discr_perc<0.01); %&monotone==1);
increment_x_new = increment_x(select);
ratio_temp_new = ratio_temp(select);
N_temp_new = length(increment_x_new);
[increment_x_ord] = sort(increment_x_new);
ratio_trimmed_mean_temp = mean(increment_x_ord(max(1,floor(N_temp_new*0.25)):floor(N_temp_new*0.75)))/increment_z
mean_ratio_temp = mean(ratio_temp_new)

display('*** FIRST DERIVATIVES ***')
select_first = (discr_perc_first<0.01); %&monotone==1);
increment_x_new_first = increment_x_first(select_first);
ratio_temp_new_first = ratio_temp_first(select_first);
N_temp_new = length(increment_x_new_first);
[increment_x_ord_first] = sort(increment_x_new_first);
ratio_trimmed_mean_temp_first = mean(increment_x_ord_first(max(1,floor(N_temp_new*0.25)):floor(N_temp_new*0.75)))/increment_z
mean_ratio_temp_first = mean(ratio_temp_new_first)


